## Replication code for Table 1 of Abrajano, Elmendorf, and Quinn.
## "Labels vs. Pictures: Treatment Mode Effects in Experiments About
##  Discrimination." Political Analysis. 
##
## Marisa Abrajano, Christopher Elmendorf, and Kevin Quinn
## 9/24/2017
##

library(lmtest)
library(sandwich)
library(msm)
library(ggplot2)

## function that does clustered SEs
## (from Hainmueller, Hopkins, & Yamomoto)
vcovCluster <- function(model, cluster){
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}


## load the long (stacked) data
load("./PhotosLong-ALL.Rda")
load("./WordsLong-ALL.Rda")

## convert choice variable into numeric
mydata.w.long$choice <- as.numeric(mydata.w.long$choice == "CandA")
mydata.p.long$choice <- as.numeric(mydata.p.long$choice == "CandA")


###########################################################################
## Hypothesis 1 (and 4 and 5)

## estimate tau (effect of text ethnicity on choice for full sample)
mydata.sub <- na.omit(mydata.w.long)
lm.tau <- lm(choice ~ candA.ethnicity, data=mydata.sub)
tab.tau <- coeftest(lm.tau, vcov=vcovCluster(lm.tau,
                                             cluster=mydata.sub$respondent.ID))


## estimate nu (effect of photo ethnicity on choice for full sample)
mydata.sub <- na.omit(mydata.p.long)
lm.nu <- lm(choice ~ candA.ethnicity, data=mydata.sub)
tab.nu <- coeftest(lm.nu, vcov=vcovCluster(lm.nu,
                                             cluster=mydata.sub$respondent.ID))



delta1.hat <- coef(lm.tau)[2] - coef(lm.nu)[2]
delta1.var <- tab.tau[2,2]^2 + tab.nu[2,2]^2
delta1.se <- sqrt(delta1.var)
delta1.z <- delta1.hat / delta1.se
delta1.pval <- 1 - pchisq(delta1.z^2, df=1)





###########################################################################
## Hypothesis 2

## estimate tau_I (effect of text ethnicity on choice for high int. mot. )
mydata.sub <- mydata.w.long[mydata.w.long$int.mot.high==TRUE,]
mydata.sub <- na.omit(mydata.sub)
lm.tauI <- lm(choice ~ candA.ethnicity, data=mydata.sub)
tab.tauI <- coeftest(lm.tauI, vcov=vcovCluster(lm.tauI,
                                             cluster=mydata.sub$respondent.ID))

## estimate tau_notI (effect of text ethnicity on choice for low int. mot. )
mydata.sub <- mydata.w.long[mydata.w.long$int.mot.high==FALSE,]
mydata.sub <- na.omit(mydata.sub)
lm.taunotI <- lm(choice ~ candA.ethnicity, data=mydata.sub)
tab.taunotI <- coeftest(lm.taunotI, vcov=vcovCluster(lm.taunotI,
                                             cluster=mydata.sub$respondent.ID))


## estimate nu_I (effect of photo ethnicity on choice for high int. mot. )
mydata.sub <- mydata.p.long[mydata.p.long$int.mot.high==TRUE,]
mydata.sub <- na.omit(mydata.sub)
lm.nuI <- lm(choice ~ candA.ethnicity, data=mydata.sub)
tab.nuI <- coeftest(lm.nuI, vcov=vcovCluster(lm.nuI,
                                             cluster=mydata.sub$respondent.ID))

## estimate nu_notI (effect of photo ethnicity on choice for low int. mot. )
mydata.sub <- mydata.p.long[mydata.p.long$int.mot.high==FALSE,]
mydata.sub <- na.omit(mydata.sub)
lm.nunotI <- lm(choice ~ candA.ethnicity, data=mydata.sub)
tab.nunotI <- coeftest(lm.nunotI, vcov=vcovCluster(lm.nunotI,
                                             cluster=mydata.sub$respondent.ID))


Idiff.hat <- coef(lm.tauI)[2] - coef(lm.nuI)[2]
Idiff.var <- tab.tauI[2,2]^2 + tab.nuI[2,2]^2
notIdiff.hat <- coef(lm.taunotI)[2] - coef(lm.nunotI)[2]
notIdiff.var <- tab.taunotI[2,2]^2 + tab.nunotI[2,2]^2
delta2.hat <- Idiff.hat - notIdiff.hat
delta2.se <- sqrt(Idiff.var + notIdiff.var)
delta2.z <- delta2.hat / delta2.se
delta2.pval <- 1 - pchisq(delta2.z^2, df=1)




###########################################################################
## Hypothesis 3

## estimate tau_E (effect of text ethnicity on choice for high ext. mot. )
mydata.sub <- mydata.w.long[mydata.w.long$int.mot.high==FALSE & mydata.w.long$ext.mot.high==TRUE,]
mydata.sub <- na.omit(mydata.sub)
lm.tauE <- lm(choice ~ candA.ethnicity, data=mydata.sub)
tab.tauE <- coeftest(lm.tauE, vcov=vcovCluster(lm.tauE,
                                             cluster=mydata.sub$respondent.ID))


## estimate nu_E (effect of photo ethnicity on choice for high ext. mot. )
mydata.sub <- mydata.p.long[mydata.p.long$int.mot.high==FALSE & mydata.p.long$ext.mot.high==TRUE,]
mydata.sub <- na.omit(mydata.sub)
lm.nuE <- lm(choice ~ candA.ethnicity, data=mydata.sub)
tab.nuE <- coeftest(lm.nuE, vcov=vcovCluster(lm.nuE,
                                             cluster=mydata.sub$respondent.ID))


delta3.hat <- coef(lm.tauE)[2] - coef(lm.nuE)[2]
delta3.var <- tab.tauE[2,2]^2 + tab.nuE[2,2]^2
delta3.se <- sqrt(delta3.var)
delta3.z <- delta3.hat / delta3.se
delta3.pval <- 1 - pchisq(delta3.z^2, df=1)








## put estimates in a nice latex table
library(xtable)

tau.hat <- coef(lm.tau)[2]
tau.se <- tab.tau[2,2]
tau.z <- tau.hat / tau.se
tau.pval <- 1 - pchisq(tau.z^2, df=1)

nu.hat <- coef(lm.nu)[2]
nu.se <- tab.nu[2,2]
nu.z <- nu.hat / nu.se
nu.pval <- 1 - pchisq(nu.z^2, df=1)


tauI.hat <- coef(lm.tauI)[2]
tauI.se <- tab.tauI[2,2]
tauI.z <- tauI.hat / tauI.se
tauI.pval <- 1 - pchisq(tauI.z^2, df=1)


nuI.hat <- coef(lm.nuI)[2]
nuI.se <- tab.nuI[2,2]
nuI.z <- nuI.hat / nuI.se
nuI.pval <- 1 - pchisq(nuI.z^2, df=1)



taunotI.hat <- coef(lm.taunotI)[2]
taunotI.se <- tab.taunotI[2,2]
taunotI.z <- taunotI.hat / taunotI.se
taunotI.pval <- 1 - pchisq(taunotI.z^2, df=1)


nunotI.hat <- coef(lm.nunotI)[2]
nunotI.se <- tab.nunotI[2,2]
nunotI.z <- nunotI.hat / nunotI.se
nunotI.pval <- 1 - pchisq(nunotI.z^2, df=1)


tauE.hat <- coef(lm.tauE)[2]
tauE.se <- tab.tauE[2,2]
tauE.z <- tauE.hat / tauE.se
tauE.pval <- 1 - pchisq(tauE.z^2, df=1)


nuE.hat <- coef(lm.nuE)[2]
nuE.se <- tab.nuE[2,2]
nuE.z <- nuE.hat / nuE.se
nuE.pval <- 1 - pchisq(nuE.z^2, df=1)

estimands <- c("tau", "nu", "delta1",
               "tauI", "nuI", "taunotI", "nunotI", "delta2",
               "tauE", "nuE", "delta3")
estimates <- c(tau.hat, nu.hat, delta1.hat,
               tauI.hat, nuI.hat, taunotI.hat, nunotI.hat, delta2.hat,
               tauE.hat, nuE.hat, delta3.hat)
ses <- c(tau.se, nu.se, delta1.se,
         tauI.se, nuI.se, taunotI.se, nunotI.se, delta2.se,
         tauE.se, nuE.se, delta3.se)
pvals <- c(tau.pval, nu.pval, delta1.pval,
               tauI.pval, nuI.pval, taunotI.pval, nunotI.pval, delta2.pval,
               tauE.pval, nuE.pval, delta3.pval)



tab.out <- data.frame(estimand=estimands, estimate=estimates,
                      se=ses, pval=pvals)

rownames(tab.out) <- NULL
print(xtable(tab.out, digits=3), include.rownames=FALSE)
